/* Integration points (quadratures)
   Jukka Toivanen, Jan Stebel
*/

#include "constpoint.h"



template<unsigned int dim>
void ConstPoint<dim>::setGSF()
{
	for (unsigned int i=0; i<dim; i++)
		GSF(i,0) = 0;
}

template<unsigned int dim>
void ConstPoint<dim>::setSF()
{
  SF(0,0) = 1;
}

template<unsigned int dim>
void ConstPoint<dim>::setX(int elem, eval_t E)
{
  x.resize(dim+1,3);
  x_a.resize(dim+1,3);
  switch (E)
  {
    case EVAL_F:
    case EVAL_FG:
      for (unsigned int i=0; i<dim+1; i++)
        for (int j=0; j<3; j++)
          x(i,j)=imesh->getCoord(IntPoint<dim>::mapDof(elem,i),j).getVal();
      break;
    case EVAL_FA:
      for (unsigned int i=0; i<dim+1; i++)
        for (int j=0; j<3; j++)
          x_a(i,j)=imesh->getCoord(IntPoint<dim>::mapDof(elem,i),j);
      break;
  }
}

template<unsigned int dim>
void ConstPoint<dim>::calcJ(eval_t E)
{
  switch (E)
  {
    case EVAL_F:
    case EVAL_FG:
      for (unsigned int i=0; i<dim; i++)
      {
        for (int j=0; j<3; j++) J(i,j) = x(1+i,j)-x(0,j);
      }
      break;
    case EVAL_FA:
      for (unsigned int i=0; i<dim; i++)
      {
        for (int j=0; j<3; j++) J_a(i,j) = x_a(1+i,j)-x_a(0,j);
      }
      break;
  }
}

template<>
inline int ConstPoint<3>::mapDof(int e, int n)
{
	return e;
}

template<>
inline int ConstPoint<3>::mapDof(int elem, int dof, int iParent)
{
	throw 1;
}


template<>
inline int ConstPoint<2>::mapDof(int e, int n)
{
	return imesh->getNvolumes() + e;
}

template<>
inline int ConstPoint<2>::mapDof(int elem, int dof, int iParent)
{
	return imesh->getFaceParentNo(elem,iParent);
}


template<>
inline int ConstPoint<1>::mapDof(int e, int n)
{
	return imesh->getNvolumes() + imesh->getNfaces() + e;
}

template<>
int ConstPoint<1>::mapDof(int elem, int dof, int iParent)
{
	return imesh->getEdgeParentNo(elem,iParent);
}


template class ConstPoint<1>;
template class ConstPoint<2>;
template class ConstPoint<3>;

